%----------------------------
% clear plots
 close all;
%---------------------------

// Christoph Gortz and John Tsoukalas
// News and Financial Intermediation in Aggregate Fluctuations, Review of Economics and Statistics

// This model is the standard two-sector setup without financial intermediaries!

// Two sector model with price and wage rigidities with permanent and transitory 
// technology shocks and separate specification of Multi Factor Productivity in I anc C sectors
// Capital is immobile across sectors and there are two investment and capital utilization decisions
// intratemporal investment adjustment cocts are incorporated a-la Huffman and Whynne
// This version does not include an Investment-specific shock v_l and a TFP transitory shock to cons sector z_l!
// Extended by risk premium shock
// Extended by correct version of Inflation Objective shock (in Taylor Rule, Price and Wage Phillips curves
// Taylor Rule similar to Smets Wouters 2005

var   labi labc lab kc ki kpi kpc mcc mci zcapi zcapc rki rkc pi lam phifi phifc c invei invec inve y pinfc pinfi w r 
z v b g qs ms sw ewma z_l v_l dy dc dinve    dpi   dw labobs robs pinfobsc pinfobsi rrate
epinfmac epinfmai spinfc spinfi
zcapif zcapcf rkcf rkif kpif kpcf kif kcf lamf phifif phifcf cf invef inveif invecf yf labf labif labcf pif wf 
spi srp
spkc spki epkc1aux epkc4aux epkc8aux epki1aux epki4aux epki8aux
ez1aux ez4aux ez8aux ev1aux ev4aux ev8aux rbc rbi qc qi tfp etfp4aux etfp8aux
;



varexo ez ev eb eg eqs em epinfc epinfi ew evl ezl espi erp
        epkc epki epkc1 epkc4 epkc8 epki1 epki4 epki8
        ez1 ez4 ez8 ev1 ev4 ev8 etfp etfp4 etfp8;
 
parameters constepinfc constepinfi constebeta cmaw cmapi cmapc alfac alfai
czcapi czcapc czcap  cdeltai cdeltac csigma chabb  
cindw cprobw cindpc cprobpc cindpi cprobpi csigl cfc
crpi crdy cry crr 
crhoz crhov crhob crhog crhoqs crhoms crhopinf crhow crhozl crhovl crhospi crhosrp
cminrho
cg cga 
cgv clandapaux clandawaux
Lss nvalueic 
crhopinfc crhopinfi
crhospkc crhospki crhotfp csadjcostc csadjcosti

// Steady state parameters
X1aux crki crkc cpi cw KcLc KiLi Lc Li css iss Kc Ki K Kibar Kcbar yss iiss icss rbcss rbiss qcss qiss
// Parameters needed to calculate ss parameters
cbeta ga 
gv clandap clandaw;


// fixed parameters
cdeltai =.025;   // depreciation rate I sector capital
cdeltac =.025;   // depreciation rate C sector capital
constebeta = 0.2607;   // equivalent to cbeta = 0.9974, constistent with 1990Q2-2011Q1 US data
constepinfc = 0.6722;  // constistent with 1990Q2-2011Q1 US data
constepinfi = 0.0245;  // constistent with 1990Q2-2011Q1 US data
alfac=.3;
alfai=.3;
cg = 0.177;            // ss nominal government spending to nominal output ratio, constistent with 1982Q4-2011Q1 US data
nvalueic = 0.399;      // average value of nominal investment to consumption ratio, constistent with 1982Q4-2011Q1 US data
Lss = 1;               // ss of total hours (by choosing appropriate value of the weight of leisure in utility)
cmapc =    0.0;          // MA price mark-up
cmapi =    0.0;          // MA price mark-up
cmaw  =   0.0;           // MA wage mark-up
cminrho = 0.0;           // intrtemporal inv. adjustment cost, zero equals crho = -1
csigma=1.0;
cfc=1.5;


// estimated parameters initialisation

cga = 0.141; 
cgv = 0.434;
csadjcost= 2.2389;	// investment adj. cost parameter
chabb=    0.6231;      // consumption habit 
cprobw=   0.6672;     // wage calvo contract
cindw=    0.1182;      // wage indexation
csigl=    0.6691;      // inverse labour supply ela
cprobpc=   0.7842;    // consumption sector price calvo probability 
cindpc=    0.0690;     // consumption sector price indexation
cprobpi=   0.7096;    // investment sector price calvo probability 
cindpi=    0.2710;     // investmnt sector price indexation
czcapi=    4.8295;    // variable I-capital utilization 
czcapc=    4.4477;    // variable C-capital utilization 
czcap =    5.75;
crpi=     1.5459;      // taylor rule inflation 
crr=      0.8423;        // taylor rule inertia  
cry=      0.0;        // taylor rule output gap 
crdy=     0.6805;        // taylor rule output gap growth
clandawaux = 1.15-1; // ss wage mark up
clandapaux = 1.15-1; // ss price mark up

crhoz=    0.7558;      // consumption sector productivity AR(1) --- growth shock
crhov=    0.1316;      // investment sector productivity AR(1) --- growth shock
crhozl=    0.0;        // consumption sector productivity AR(1) --- stationary shock
crhovl=    0.0;        // investment sector productivity AR(1) --- stationary shock
crhob=    0.9177;      // preference shock AR(1)
crhog=    0.9921;      // gov spending shock AR(1)
crhoqs=   0.0;         // installation shock AR(1)
crhoms=   0;            // monetary policy shock AR(1)
crhopinfi= 0.8968;     // price mark-up shock AR(1) in I sector
crhopinfc= 0.0456;     // price mark-up shock AR(1) in C sector
crhow=    0.0384;      // wage mark-up shock AR(1)
crhospkc = 0.8514;     // quality of physical capital shock in C sector
crhospki = 0.0682;     // quality of physical capital shock in I sector
crhospi = 0.0;         // persistence of inflation objective shock
crhopinf = 0.0;
crhozl  = 0.0;
crhovl = 0.0;
crhosrp = 0.0;
crhotfp = 0.95;      // Persistence common TFP


// Steady state parameters
// Initial values according to steady state with rho = -1
// parameters needed for this section:
cbeta = 1/(1+constebeta/100);
ga = cga/100 ;
gv = cgv/100 ;
clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;




model (linear); 

// Auxiliary parameter
//# crho = cminrho;  
# crho = 1/(cminrho - 1);  // rho = -rho
//# cga = constepinfi - constepinfc -  (alfac-1)/(1-alfai)*cgv;
//# ga = cga/100;


// flexible economy

          
          invef = clandap*(alfai*kif + (1-alfai)*labif + v_l);
          invef = (iiss^-crho + icss^-crho)^-1*iiss^-crho*inveif + (iiss^-crho + icss^-crho)^-1*icss^-crho*invecf;
          cf + crki*Kibar/css*exp(-(1/(1-alfai))*gv)*zcapif + crkc*Kcbar/css*exp(-(1/(1-alfai))*gv)*zcapcf = clandap*(alfac*kcf + (1-alfac)*labcf + z_l);
	      z_l - alfac*labcf + alfac*kcf = wf;
          z_l + (1 - alfac)*labcf + (alfac-1)*kcf = rkcf;
          v_l - alfai*labif + alfai*kif + pif = wf;
          v_l + (1 - alfai)*labif + (alfai-1)*kif + pif = rkif;
          
      lamf = (chabb*cbeta*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf(+1)
        - ((exp(2*(ga+(alfac/(1-alfai))*gv)) + chabb^2*cbeta)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf
        + (chabb*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*cf(-1)
        + ( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhoz - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*z
	    + b
        +  (alfac/(1-alfai))*( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhov - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*v;

        // normalization: ( exp(ga+(alfac/(1-alfai))*gv) - chabb*cbeta*crhob )/(exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*b
        
        phifif = (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv)*( phifif(+1) - z(+1) - 1/(1-alfai)*v(+1) + spki(+1))
        + (1 - (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv))*( lamf(+1) - z(+1) - 1/(1-alfai)*v(+1) + rkif(+1) + zcapif(+1) + spki(+1) );

        phifcf = (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv)*( phifcf(+1) - z(+1) - 1/(1-alfai)*v(+1)+ spkc(+1) )
        + (1 - (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv))*( lamf(+1) - z(+1) - 1/(1-alfai)*v(+1) + rkcf(+1) + zcapcf(+1) + spkc(+1) );
       
        rkif = czcapi*zcapif;

        rkcf = czcapc*zcapcf;
        
       
        lamf = -pif + phifif + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(inveif - inveif(-1) + z + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( inveif(+1) - inveif + z(+1) +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1*(icss^-crho*invecf+iiss^-crho*inveif)-inveif);



        lamf = -pif + phifcf + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invecf - invecf(-1) + z + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invecf(+1) - invecf + z(+1) +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1*(icss^-crho*invecf+iiss^-crho*inveif)-invecf);
        
        
	    kif =  kpif(-1)+ zcapif + spki - 1/(1-alfai)*v;

        kcf =  kpcf(-1)+ zcapcf + spkc - 1/(1-alfai)*v;
        
        kpif =  (1-cdeltai)*exp(-(1/(1-alfai))*gv)*( kpif(-1) + spki - (1/(1-alfai))*v ) + ( 1 - (1-cdeltai)*exp(-(1/(1-alfai))*gv))*(qs + inveif);
        
        kpcf =  (1-cdeltac)*exp(-(1/(1-alfai))*gv)*( kpcf(-1) + spkc - (1/(1-alfai))*v ) + ( 1 - (1-cdeltac)*exp(-(1/(1-alfai))*gv))*(qs + invecf);

        wf =  b + csigl*labf - lamf;
        
        (1-cg)*yf = css/yss*cf+cpi*iss/yss*(invef + pif) + (1-cg)*g ;
        
        labf = Lc/Lss*labcf + Li/Lss*labif;
        
     
       

// sticky price - wage economy


         

          // consumption and investment sector output and cost minimization
          
          inve = clandap*(alfai*ki + (1-alfai)*labi + v_l);
          inve = (iiss^-crho + icss^-crho)^-1*iiss^-crho*invei + (iiss^-crho + icss^-crho)^-1*icss^-crho*invec;
          //inve = (iiss/iss)*invei + (icss/iss)*invec;
          c  + crki*Kibar/css*exp(-(1/(1-alfai))*gv)*zcapi + crkc*Kcbar/css*exp(-(1/(1-alfai))*gv)*zcapc = clandap*(alfac*kc + (1-alfac)*labc + z_l);

	      mcc =  alfac*rkc+(1-alfac)*w - z_l;
          mci =  alfai*rki+(1-alfai)*w - pi - v_l;
          rkc - w = labc - kc;
          rki - w = labi - ki;


               
          
         // consumption and investment sector optimal pricing --- Philips curves
          
           pinfc - spi =  (cbeta/(1+cindpc*cbeta))*(pinfc(+1)-spi) + (cindpc/(1+cindpc*cbeta))*(pinfc(-1)-spi) 
               +( (1-cprobpc)*(1-cbeta*cprobpc)/((1+cindpc*cbeta)*cprobpc) )* mcc  +  spinfc ; 
          
               
           pinfi - spi =  (cbeta/(1+cindpi*cbeta))*(pinfi(+1)-spi) + (cindpi/(1+cindpi*cbeta))*(pinfi(-1)-spi) 
               +( (1-cprobpi)*(1-cbeta*cprobpi)/((1+cindpi*cbeta)*cprobpi) )* mci  + spinfi ; 
          
          // normalization of shocks:
          // spinf* = ( (1-cprobpi)*(1-cbeta*cprobpi)/((1+cindpi*cbeta)*cprobpi) )* spinf ; 
            
        // Households
        
        lam = (chabb*cbeta*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c(+1)
        - ((exp(2*(ga+(alfac/(1-alfai))*gv)) + chabb^2*cbeta)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c
        + (chabb*exp(ga+(alfac/(1-alfai))*gv)/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb)))*c(-1)
        + ( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhoz - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*z
	    + b
        +  (alfac/(1-alfai))*( exp(ga+(alfac/(1-alfai))*gv)*chabb*cbeta*crhov - chabb*exp(ga+(alfac/(1-alfai))*gv) )/((exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*(exp(ga+(alfac/(1-alfai))*gv)-chabb))*v;
        
	    //normalization of shocks:
        //b* =  ( exp(ga+(alfac/(1-alfai))*gv) - chabb*cbeta*crhob )/(exp(ga+(alfac/(1-alfai))*gv)-chabb*cbeta)*b

       
        lam = r + srp + lam(+1) - z(+1) - v(+1)*alfac/(1-alfai) - pinfc(+1);
        
        phifi = (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv)*( phifi(+1) - 1/(1-alfai)*v(+1) +spki(+1) )
        + (1 - (1-cdeltai)*cbeta*exp(-(1/(1-alfai))*gv))*( lam(+1) - 1/(1-alfai)*v(+1) + rki(+1) + zcapi(+1) + spki(+1) );

        phifc = (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv)*( phifc(+1) - 1/(1-alfai)*v(+1) +spkc(+1) )
        + (1 - (1-cdeltac)*cbeta*exp(-(1/(1-alfai))*gv))*( lam(+1) - 1/(1-alfai)*v(+1) + rkc(+1) + zcapc(+1) + spkc(+1) );


        
        rki = czcapi*zcapi;
        rkc = czcapc*zcapc;

        //zcapi =  (1/(czcap/(1-czcap)))* rki ;
        //zcapc =  (1/(czcap/(1-czcap)))* rkc ;

         
        lam = -pi + phifi + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invei - invei(-1) + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invei(+1) - invei +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1* (icss^-crho*invec+iiss^-crho*invei)-invei);

        lam = -pi + phifc + qs - exp(2*((1/(1-alfai))*gv))*csadjcostc*(invec - invec(-1) + 1/(1-alfai)*v) 
        + cbeta*exp(2*((1/(1-alfai))*gv))*csadjcostc*( invec(+1) - invec +  1/(1-alfai)*v(+1) ) -(1+crho)*((iiss^-crho + icss^-crho)^-1* (icss^-crho*invec+iiss^-crho*invei)-invec);
        
 
        
	    ki =  kpi(-1)+ zcapi + spki - 1/(1-alfai)*v;

        kc =  kpc(-1)+ zcapc + spkc - 1/(1-alfai)*v;
        
        kpi =  (1-cdeltai)*exp(-(1/(1-alfai))*gv)*( kpi(-1) + spki - (1/(1-alfai))*v ) + ( 1 - (1-cdeltai)*exp(-(1/(1-alfai))*gv))*(qs + invei);

        kpc =  (1-cdeltac)*exp(-(1/(1-alfai))*gv)*( kpc(-1) + spkc - (1/(1-alfai))*v ) + ( 1 - (1-cdeltac)*exp(-(1/(1-alfai))*gv))*(qs + invec);
             
	       
	      w =  (1/(1+cbeta))*w(-1)
               +(cbeta/(1+cbeta))*w(+1)
               - ( (1-cprobw*cbeta)*(1-cprobw)/(cprobw*(1+cbeta)*(1+csigl*(1+1/(clandawaux)))))* (w - csigl*lab - b + lam)
               + (cindw/(1+cbeta))*(pinfc(-1)-spi)
               - ((1+cbeta*cindw)/(1+cbeta))*(pinfc-spi)
               + (cbeta/(1+cbeta))*(pinfc(+1)-spi)
               + (cindw/(1+cbeta))*( z(-1) + alfac/(1-alfai)*v(-1) )
               - (1+cbeta*cindw-crhoz*cbeta)/(1+cbeta)*z
               - ((1+cbeta*cindw-crhov*cbeta)/(1+cbeta))*(alfac/(1-alfai))*v
               + sw ;
               
      // normalization of shocks:
      //       sw* = ( (1-cprobw*cbeta)*(1-cprobw)/(cprobw*(1+cbeta)*(1+csigl*(1+1/(clandaw-1)))))* sw 
               
      lab = Lc/Lss*labc + Li/Lss*labi;

               
            y = (css/yss)*c+(cpi*iss/yss)*(inve+ pi) +g ;

               

     
      // Taylor rule    
   
  
	     r =  crpi*(1-crr)*pinfc
               +cry*(1-crr)*0
               +crdy*(1-crr)*(y-y(-1))
               +crr*r(-1)
              +ms  ;

        

    pinfi - pinfc  = pi - pi(-1); 

        //
    	rbc = 1/(crkc+qcss*(1-cdeltac)) * ( crkc*(rkc+zcapc) + qcss*(1-cdeltac)*(qc) ) -qc(-1) +spkc + z - (1-alfac)/(1-alfai)*v;

        // Stochastic return on assets for bank in I sector
    	rbi = 1/(crki+qiss*(1-cdeltai)) * ( crki*(rki+zcapi) + qiss*(1-cdeltai)*(qi) ) -qi(-1) +spki + z - (1-alfac)/(1-alfai)*v;

        // Price of capital in C sector
		qc = exp(2*1/(1-alfai)*gv)*csadjcostc * ( invec-invec(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (invec(+1)-invec+1/(1-alfai)*v(+1)) - qs + pi + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invec + iiss^(-crho)*invei) - invec );

        // Price of capital in I sector
		qi = exp(2*1/(1-alfai)*gv)*csadjcostc * ( invei-invei(-1)+1/(1-alfai)*v) - cbeta*exp(2*1/(1-alfai)*gv)*csadjcostc * (invei(+1)-invei+1/(1-alfai)*v(+1)) - qs + pi + (1+crho) * ( (iiss^(-crho)+icss^(-crho))^(-1) * (icss^(-crho)*invec + iiss^(-crho)*invei) - invei );



  //dpi + (constepinfc - constepinfi) = pinfi  - pinfc;
// this implies the restriction
// constepinfi - constepinfc = (cga+ (alfac-1)/(1-alfai)*cgv)
// so initial parameter values need to respect this
// in order to solve the SS 
      


    // TFP growth shock to cons sector

        z = crhoz*z(-1)  + ez  + ez1aux(-1) + ez4aux(-4) + ez8aux(-8);
        ez1aux = ez1;
        ez4aux = ez4;
        ez8aux = ez8;
    
    // TFP stationary shock to cons sector
      
        z_l = crhozl*z_l(-1)  + ezl;
         
   // Investment-specific sector growth shock

        v = crhov*v(-1)  + ev + ev1aux(-1) + ev4aux(-4) + ev8aux(-8);
        ev1aux = ev1;
        ev4aux = ev4;
        ev8aux = ev8;

    // Investment-specific stationary shock 
      
         v_l = crhovl*v_l(-1)  + evl;

   
   //  Preference, government, installation
   //  monetary policy, price markup, wage markup
	           
        b = crhob*b(-1) + eb;
	      
        g = crhog*g(-1) + eg;
        
        qs = crhoqs*qs(-1) + eqs;
	     
        ms = crhoms*ms(-1) + em;
	    
        spinfc = crhopinfc*spinfc(-1) + epinfmac - cmapc*epinfmac(-1);
	          epinfmac=epinfc;

        spinfi = crhopinfi*spinfi(-1) + epinfmai - cmapi*epinfmai(-1);
	          epinfmai=epinfi;
	    
        sw = crhow*sw(-1) + ewma - cmaw*ewma(-1) ;
	          ewma=ew; 
    // Shock to inflation objective
	      spi = crhospi*spi(-1) + espi;
    // Risk premium shock
        srp = crhosrp*srp(-1) + erp;

        // Shock to the quality of physical capital in C and I sector
        spkc = crhospkc*spkc(-1) + epkc + epkc1aux(-1) + epkc4aux(-4) + epkc8aux(-8);
        epkc1aux = epkc1;
        //epkc2aux = epkc2;
        //epkc3aux = epkc3;
        epkc4aux = epkc4;
        epkc8aux = epkc8;
        spki = crhospki*spki(-1) + epki + epki1aux(-1) + epki4aux(-4) + epki8aux(-8);
        epki1aux = epki1;
        //epki2aux = epki2;
        //epki3aux = epki3;
        epki4aux = epki4;
        epki8aux = epki8;

// TFP shock common to both sectors
        tfp = crhotfp*tfp(-1) + etfp + etfp4aux(-4) + etfp8aux(-8);
        etfp4aux = etfp4;
        etfp8aux = etfp8;



// measurement equations
/*
dy=y-y(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv) ;
dc=c-c(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv)  ;
dinve=inve-inve(-1) + 1/(1-alfai)*v + ( 1/(1-alfai)*cgv) ;
dpi = pi - pi(-1) + z + (alfac-1)/(1-alfai)*v + (cga+ (alfac-1)/(1-alfai)*cgv) ;
dw=w-w(-1) + z + alfac/(1-alfai)*v + (cga + alfac/(1-alfai)*cgv)  ;
pinfobsc = 1*(pinfc) + constepinfc ;
pinfobsi = 1*(pinfi) + constepinfi ;
robs =    r + log(constepinfc/100 + 1) - log(cbeta) + (cga + alfac/(1-alfai)*cgv) ;
labobs = lab + log(Lss) ;
rrate = r - pinfc(+1);
*/



// NEW measurement equations (using demeaned data)

dy=y-y(-1) + z + alfac/(1-alfai)*v  ;
dc=c-c(-1) + z + alfac/(1-alfai)*v   ;
dinve=inve-inve(-1) + 1/(1-alfai)*v  ;
dpi = pi - pi(-1) + z + (alfac-1)/(1-alfai)*v  ;
dw=w-w(-1) + z + alfac/(1-alfai)*v   ;
pinfobsc = 1*(pinfc)  ;
pinfobsi = 1*(pinfi)  ;
robs =    r  ;
labobs = lab + log(Lss) ;
rrate = r - pinfc(+1);



end; 

//steady;

shocks;

var ez;
stderr 0.1611;
var ev;
stderr 1.8366;
var eb;
stderr 1.3419;
var eg;
stderr 0.5019;
var eqs;
stderr 0;
var em;
stderr 0.1175;
var epinfc;
stderr 0.5929;
var epinfi;
stderr 0.2150;
var ew;
stderr 0.3618;
var espi;
stderr 0;
var erp;
stderr 0;
var ezl;
stderr 0.0;
var evl;
stderr 0.0;
var epki;
stderr 2.3752;
var epkc;
stderr 0.2835;
var epkc1;
stderr 0; 
var epkc4; 
stderr 0;
var epkc8; 
stderr 0;  
var epki1; 
stderr 0.0;
var epki4;
stderr 0.0;
var epki8;
stderr 0.0;
var ez1;
stderr 0;
var ez4;
stderr 0.1095;
var ez8;
stderr 0.1917;
var ev1;
stderr 0.0;
var ev4;
stderr 0.2352;
var ev8;
stderr 0.7066;


var etfp;
stderr 0.0;
var etfp4;
stderr 0.0/sqrt(2);
var etfp8;
stderr 0.0/sqrt(2);

end;



estimated_params;

// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF


stderr ez,0.132,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr ev,1.2186,0.001,100,INV_GAMMA_PDF,0.5,2.0;
stderr eb,0.247,0.025,150,INV_GAMMA_PDF,0.1,2.0;
stderr eg,0.39,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr eqs,1.27,0.01,120,INV_GAMMA_PDF,0.1,2.0;
stderr em,0.113,0.01,80,INV_GAMMA_PDF,0.1,2.0;
//stderr espi,0.20,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr epinfc,0.57,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr epinfi,0.27,0.01,80,INV_GAMMA_PDF,0.1,2.0;
stderr ew,0.42,0.01,20,INV_GAMMA_PDF,0.1,2.0;
//stderr enc,0.260,0.01,100,INV_GAMMA_PDF,0.1,2.0;
//stderr eni,0.19,0.01,100,INV_GAMMA_PDF,0.1,2.0;
stderr epkc,0.860,0.01,100,INV_GAMMA_PDF,0.5,2.0;
stderr epki,1.15370,0.01,100,INV_GAMMA_PDF,0.5,2.0;
//stderr ethetab,0.001,0.0001,10,INV_GAMMA_PDF,0.001,2.0;
//stderr ezl,0.57,0.01,40,INV_GAMMA_PDF,0.5,2.0;
//stderr evl,0.95,0.01,40,INV_GAMMA_PDF,0.5,2.0;
//stderr esp1, 0.05,INV_GAMMA_PDF,0.05,0.01;
//stderr esp2, 0.05,INV_GAMMA_PDF,0.05,0.01;
//stderr epkc1,0.25,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epkc2,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epkc3,0.15,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epkc4,0.6,0.01,100,INV_GAMMA_PDF,0.1/sqrt(2),2.0;
//stderr epkc8,0.95,0.01,100,INV_GAMMA_PDF,0.1/sqrt(2),2.0;
//stderr epki1,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epki2,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epki3,0.05,0.01,100,INV_GAMMA_PDF,0.1/sqrt(5),2.0;
//stderr epki4,1.03,0.01,100,INV_GAMMA_PDF,0.1/sqrt(2),2.0;
//stderr epki8,0.85,0.01,100,INV_GAMMA_PDF,0.1/sqrt(2),2.0;

//stderr etfp,0.139,0.01,100,INV_GAMMA_PDF,0.5,2.0;



//stderr ez1,0.05,0.01,100,INV_GAMMA_PDF,0.5/sqrt(1),2.0;
stderr ez4,0.315,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
stderr ez8,0.1675,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr ev1,0.05,0.01,100,INV_GAMMA_PDF,0.5/sqrt(1),2.0;
stderr ev4,0.3185,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
stderr ev8,0.45,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;

//stderr etfp4,0.3285,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;
//stderr etfp8,0.3185,0.01,100,INV_GAMMA_PDF,0.5/sqrt(2),2.0;


crhoz,0.41,.01,.9999,BETA_PDF,0.4,0.20;
crhov,.6 ,.01,.9999,BETA_PDF,0.4,0.20;
crhob,.949,.01,.9999,BETA_PDF,0.6,0.20;
crhog,.9,.01,.9999,BETA_PDF,0.6,0.20;
//crhoqs,.7,.01,.9999,BETA_PDF,0.6,0.20;
//crhoms,.38,.0000000001,.9999,BETA_PDF,0.4,0.20;
//crhospi,.95,.0000000001,.9999,BETA_PDF,0.4,0.20;
crhopinfc,.637,.01,.9999,BETA_PDF,0.6,0.20;
crhopinfi,.92,.01,.9999,BETA_PDF,0.6,0.20;
crhow,.84,.001,.9999,BETA_PDF,0.6,0.20;
//crhosnc,0.82 ,.01,.9999,BETA_PDF,0.5,0.20;
//crhosni,0.55 ,.01,.9999,BETA_PDF,0.5,0.20;
crhospkc,0.8 ,.01,.9999,BETA_PDF,0.6,0.20;
crhospki,0.9 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhozl,.936 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhovl,.532 ,.01,.9999,BETA_PDF,0.6,0.20;
//crhothetab,0.7 ,.01,.9999,BETA_PDF,0.5,0.20;
//crhotfp,0.837 ,.01,.9999,BETA_PDF,0.6,0.20;


//cmapc,.428,0.01,.9999,BETA_PDF,0.5,0.2;
//cmapi,.8190,0.01,.9999,BETA_PDF,0.5,0.2;
//cmaw,.93,0.01,.9999,BETA_PDF,0.5,0.2;


csadjcostc,1.4,0.2,50,GAMMA_PDF,4,1.0;
//csadjcosti,1.4,0.2,50,GAMMA_PDF,4,1.0;

chabb,0.537,0.001,0.99,BETA_PDF,0.5,0.1;
cprobw,0.69,0.01,0.99,BETA_PDF,0.66,0.1;
csigl,1.9317,0.025,10,GAMMA_PDF,2,0.75;
cprobpc,0.681,0.01,0.99,BETA_PDF,0.66,0.10;
cprobpi,0.7181,0.01,0.99,BETA_PDF,0.66,0.10;
cindw,0.20,0.000000001,0.99,BETA_PDF,0.5,0.15;
cindpc,0.1683,0.000000001,0.99,BETA_PDF,0.5,0.15;
cindpi,0.5429,0.000000001,0.99,BETA_PDF,0.5,0.15;
czcapi,4.68382,0.000001,50,GAMMA_PDF,5.0,1.0;
czcapc,5.447855,0.000001,50,GAMMA_PDF,5.0,1.0;
//czcap,12.2,0.000001,50,GAMMA_PDF,5.0,1.0;
//czcap,0.8,0.000001,1,BETA_PDF,0.5,0.1;

crpi,1.4457504,1.0,5,NORMAL_PDF,1.7,0.30;
crr,0.83,0.1,0.99,BETA_PDF,0.60,0.20;
//cry,0.09,0.0000000001,0.8,NORMAL_PDF,0.125,0.05;
crdy,0.69,0.000001,0.9,NORMAL_PDF,0.25,0.1;
//cry,0.0467,0.0000000001,0.8,NORMAL_PDF,0.25,0.1;

//cminrho,0.88, 0.00000000001,0.999999,BETA_PDF,0.5,0.20;


end;

varobs dc dinve dy labobs dw  pinfobsc pinfobsi robs ; //dequity spread dcredit;
// For Estimation
estimation(datafile= USdata_10, mode_compute=0, mode_file = twosb_v10a_NoFI_mode,  first_obs=1, presample=1,  prior_trunc=1e-50, lik_init=1, prefilter=0, mh_replic=0,   mh_nblocks=2,mh_jscale=0.25,mh_drop=0.2,  nodiagnostic ) dc;
//estimation(datafile= USdata_6, mode_compute=4, mode_file = twosb2v6_mode, mode_check, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=0, mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2, nodiagnostic) dc;
//estimation(datafile= USdata_6, mode_compute=0, mode_file = twosb2v6_mode, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=50000, mh_nblocks=2,mh_jscale=0.26,mh_drop=0.2) dc;
//load_mh_file,

// MDD after 50000:
// MDD 

//L0gL = -471.14
//Poster = -560

// For IRFs
//estimation(datafile= USdata_6, mode_compute=0, mode_file = twosb2v6_mode, prior_trunc=1e-50, first_obs=1, presample=1, lik_init=1, prefilter=0, mh_replic=0, mh_nblocks=2,mh_jscale=0.25,mh_drop=0.2) dc;
//stoch_simul(irf = 0) dy dc dinve dw pinfobsc pinfobsi robs labobs spreadobsi spreadobsc;

/*
//mode_file = twosb_v10a_NoFI_mode
cprobw=   0.01; 
cindw = 0.01;
cprobpc=   0.01; 
//cprobpi = 0.01;
clandawaux = 1.01-1; // ss wage mark up
clandapaux = 1.01-1; // ss wage mark up


clandap = 1 + clandapaux;
clandaw = 1 + clandawaux;

cindpc = 0.01;
//cindpi = 0.01;

//csadjcostc = 0.01;

*/
stoch_simul( irf = 0)   dc dinve dy labobs dw  pinfobsc pinfobsi robs rkc rki;//dc dinve dy labobs labi labc invec invei dw  pinfobsc pinfobsi robs rbc rbi rrate qc qi  ;//y inve c labobs labi labc spreadobsc spreadobsi robs pinfobsc dpi nc ni;